Phase-field model with reduced interface diffuseness 
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CZ3 i We minimized the interface diffuseness in the phase-field models by introducing the parabolic 

double-well potential and localizing the solute redistribution (or latent heat release) into a narrow 
region within the phase-field interface. In spite of the parabolic potential with cusps, highly localized 
solute redistribution and discontinuous diffusivity function adopted in this model, it works remark- 
ably well in numerical computations. The computations on dendritic solidification of an one-sided 
system yield quantitatively the same results with the anti-trapping model [A. Karma, Phys. Rev. 
Lett. 87, 115701 (2001)], indicating the anomalous interfacial effects can be effectively suppressed. 
This approach can be easily extended to the multi-components or multi-phases system. 
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I. INTRODUCTION 



Solidification processes in materials produce complex microstructural patterns in the interface morphology and 
composition distribution profile. The prediction and control of solidification microstructures are important because 
the properties of the materials are usually dictated by the microstructure. However the computational prediction of 
the microstructure has been one of the formidable free boundary problems. Traditional front tracking methods solving 
the equations of motion of the sharp interface and keeping tracking the interface position in every time steps have been 
used, but appeared to be not suitable especially when the geometry is complex or in 3D. The phase-field model (PFM) 
has emerged as an effective tool in describing the complex pattern evolution Q. In this approach, the phase field <f> 
is introduced, varying from one value in solid to another value in the liquid across a spatially diffuse interface with a 
width 2£. The equations governing the phase field and diffusion field in PFM follow naturally from the definition of 
I ■ the free energy functional of the system and requiring it to decrease monotonically in time. In particular, Karma and 
Rappels' asymptotic analysis M at the thin interface limit (2£ <C D/V, D: thermal diffusivity, V: interface velocity) 
improved greatly the computational efficiency and made the PFM to be a quantitative computational tool under the 
real experimental conditions, even though their approach was restricted to the case of pure materials' solidification 
with equal thermal diffusivity and heat capacity in solid and liquid phases. 

Several PFMs for alloys |3[-|i| also have been proposed. To use the alloy PFMs as the computational tools that can 
yield quantitative results under realistic conditions, it is highly desirable to determine the parameters in the PFMs 
at the thin interface limit, not at the sharp interface limit. As shown by Almgren ||, however, Karma and Rappel's 
asymptotic analysis can not be straightforwardly extended to the alloy PFMs because the unequal solute diffusivitics 
and solute partitioning between solid and liquid give rise to anomalous interface effects; chemical potential jump, 
surface stretching and surface diffusion. If the heat capacities and thermal diffusivities are different in solid and 
liquid, the PFMs for pure materials suffer from the same anomalous interface effects. Recently, Karma Q proposed a 
new PFM for alloy solidification, which is based on i) introduction of an anti-trapping flux in the diffusion equation 
and ii) assumption of negligible solute diffusivity in solid (Dg), compared with that in liquid (-Dl). The anti-trapping 
flux could be manipulated to not only cancel out the solute trapping current across the interface at low interface 
velocity, but also eliminate all the anomalous interface effects simultaneously, together with proper choices for some 
auxiliarly functions. Furthermore the assumption D$ <C Dl, as in the our previous approach of 1-D thin interface 
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analysis make the steady diffusion equation to be integrated without unknown constant, which enable us to find 
the relationship between the phase-field mobility and the kinetic coefficient at the thin interface limit. 

In this study we propose a simple method for getting rid of the anomalous interface effects in the PFMs, which 
can be straightforwardly extended to the multi-phases or multi-components system. The anomalous interface effects 
in PFMs originate from the diffuseness of the interface in the diffusion equation, irrespective of the diffuseness in the 
phase-field equation; they all vanish when the sharp interface limit is taken for the diffusion equation. Furthermore, 
the interface width in PFMs can be defined differently each other for diffusion field and for phase field; 2£d and 
2£ p , respectively. Our approach thus is based on the decoupling the interface widths for both fields within the thin 
interface limit and pushing the interface width for diffusion field to be minimized. Mathematically this corresponds 
to the condition Dl/V ^> 2£ p ^> 2^ — > 0, in contrast with the previous thin interface limit Dl/V ^> 2£^ ~ 2£ p 
. In spirit, present approach is somewhat similar to the pseudo-front tracking method Q proposed recently, where 
the diffusion equation is solved on the sharp interface and the curvature for the Gibbs-Thomson equation is obtained 
from the diffuse interface. In followings, we describe our approach in detail and then present computation results on 
the dendritic growth. 



II. MODEL 



We start from finding the governing equations unified for both the pure thermal case and the pure solutal case. 
With definitions </> = — 1 at liquid, <fi = 1 at solid and — 1 < (f> < 1 at the interface region, the phase-field equation 0] 
is given by 

Td^ = W 2 V 2 ^-d ( ,[f^) + Xug(^)}, (1) 

where A = aiW/do, do is the capillary length, /(</>) is the double well potential with the minima at (j> = ±1, g(<j>) 
is an odd function of <j> with </(±l) = 0, u is the dimensionless field defined as u — (T — T m ) / (AH m /Cp) for pure 
materials and u = (cl — c- L )/{c L — c e s ) for alloys, where T m is the melting point, AH m is the latent heat, Cjf and 
Cp are the specific heats of solid and liquid, respectively and c e L and c e s are the equilibrium composition of solid and 
liquid at a given temperature, respectively. When the constant a\ is taken as 
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V2/(0o)#o. (2) 



Eq. (|l|) recover the Gibbs-Thomson equation u = —d$k at equilibrium state, where k is the curvature of the interface. 
For solidification of pure materials, the diffusion equation can be written as 



, S l + H<f>) , „ L l-h{<f,), 



d t H = D L V ■ q{4>)[C b p f^+C% ^]VT, ( 3 ) 

where the local enthalpy density H is given by H = Hs[l + h(<f))]/2 + Hl[1 — h(<j>)]/2, H$ and Hl are the enthalpy 
densities of solid and liquid as functions of temperature, respectively, h(<f>) is an odd function of h(±l) = ±1, C p is 
the specific heat given by C p = Cr [1 + h(<p)]/2 + Cp [1 — h(4>)]/2, the diffusivity function q((f>) is defined to satisfy 
<?(— 1) = 1 and q(l) — Ds/Dl, where D$ and Dl are the thermal diffusivities in solid and liquid, respectively. For 
isothermal solidification of binary alloys, the solute diffusion equation [|],^| can be written as 

d t c = D L V • ^(HMv, + ^I^Vc,], (4) 

where the local composition (mole fraction) c is given by c = cg[l+/i(0)]/2+ci,[l — h((f>)]/2, Dl is the solute diffusivity 
in liquid and the corresponding q((f>) of solute diffusivity to that of the thermal diffusivity. Although the chemical 
potential that can be defined everywhere may be used instead of cs and cl as in it is more convenient to use 

Eq. (||) directly, as long as cl and cs are defined at the whole space of the system. With Hs = H™ + Cg(T - T m ), 
Hl = H™ + Cp(T — T m ) and the usual condition cs/cl = c e s /c e L [@J|], the diffusion equations (||) and @ can be 
written as an unified form; 

^M(0) - ^h(cf>)] = D L V ■ q(4>)A(0)Vu, (5) 

where A(<f>) = [l + k-(l-k)h(4>)]/2, and k = C^/C^ for solidification of pure materials and k — c e s j c% for isothermal 
solidification of binary alloys. Note that in fact the final forms of the governing equations (111) and (H) in this study 
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are fundamentally identical with those in the previous studies except the anti-trapping term in the diffusion 

equation of Ref. . 

The interface width 2£ p of phase field is determined by the form of the double well potential f(<f>) in Eq. ([!]). In 
this study a parabolic potential f(<p) = |1 — <fi 2 \/2 is adopted, instead of the traditional form (1 — (f> 2 ) 2 . The parabolic 
potential is not new one, but has been used in PFMs HJlIJ. Comparing with the fourth order potential, the parabolic 
potential gives several benefits as follows; i) The equilibrium phase- field profile is given by 4>o — ~ sva.(x/W). Outside 
of the interface width 2£ p = ttW , the phase state becomes completely either solid or liquid, which is contrast to long 
smearing of phase field into solid and liquid phases in the PFMs with fourth order potential. This clear cut of interface 
region suppresses the spurious interaction between interfaces during computation on the system with complex interface 
geometry; For example, when two solid/liquid interfaces are approaching each other, the two interfaces remains intact 
until the distance between two interface becomes less than irW. However, in case of the fourth order potential, 
two approaching interfaces interact each other with spurious long-range attraction and then merges easily into an 
interface, due to the overlapping of smearing phase fields at liquid between two solids, ii) Below we will introduce a 
highly localized h{<j)) form to minimize the anomalous interface effects. However coupling the h(4>) with the fourth 
order potential makes the phase-field profile to be less stable in computation with increased grid size, compared with 
the coupling with the parabolic potential, iii) Importantly, the parabolic potential is more convenient in extending 
the phase-field model to multiphase-fields case than the fourth order potential, because combination of the parabolic 
potentials gives local maxima at the triple junction p0| . Indeed the multi-phases-field model with the parabolic 
potentials is extensively used recently for quantitative computation of the microstructure under the real experimental 
conditions H]. 

The interface width 2£^ for diffusion field is determined by the function h(cj)) in Eq. (||). Here we adopt a simple 
form; 

r -1, for 4> < -(/>*; 
h(4>) = \ <j>/ct>*, for -</,*<(/,< +0*; (6) 
l+l, for 4>>+<t>*. 

This form of h(4>) represents that solute redistribution or latent heat release takes place in a localized region — 0* < 
(f> < +cf>* within the interfacial region (—1 < <j) < 1) of the phase field. The <f>* represents the degree of localization, 
the smaller the value the more localization of solute redistribution. The <j>* value defines a new interface width 2£d 
with the relationship 2£d = 2W smT 1 4>* , and the width is independent of the interface width 2£„ of phase field. If wc 
take <fi* = 1, then it follows h{4>) = 0, which has been widely used in previous approaches ^,||- All the anomalous 
interface effects become smaller with decreasing 2^ and completely disappear at the limit 2^ — > 0. In reality of 
numerical computation, however, vanishing 2£<j can not be taken as long as a finite mesh size is used in computations: 
Consider ID interface moving steadily in computation with mesh size Ax. If we set 2^ (over where <j> changes from 
— cf>* to +(/>*) to be smaller than Ax, there can be the moment when the phase-field values on all the grid points are 
outside the range —(f)* < <j) < +<f>*. During this moment, the system loses the existence of the interface region for 
diffusion field. Even though this phenomena does not lead to the violation of the mass or energy conservation, it plays 
the role of a very large noise source in computations. Thus the best choice of the width 2^ appears to be 2^ = Ax. 

As long as we can not take the limit 2^ — > 0, we should take care of the anomalous interface effects. The conditions 
for vanishing anomalous interface effects |6|,f7j are given by 

f+id 

[l-h(<fa)]dn= I [l + h(<p )]dx, (7) 

J -id 

for the interface stretching, 



-id Jo 



(8) 

J q{<P )A{(t) y 



'-id q(<h)A(<h) 

for the diffusion potential (chemical potential, concentration or temperature) jump and 



q(cjy )A(4> )dx= / [1 - q{c/> )A(<f> )}dx, (9) 

-id Jo 

for the interface diffusion. </>o is the equilibrium phase-field profile given by 0o = — sin(x/W) for the parabolic 
potential in this study. It should be noted that for the last condition (^) we assumed q{<f>o) — Ds/D^ — > at 
x < —£d- The first condition (Q) is satisfied with the form (^) for h(<po)- The second condition also can be satisfied by 
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simply putting q((f>o)A(<f>o) = 1 at x > — £<j. The last condition however is violated, as long as the conditions (f7j) and 
(||) are satisfied. How large the anomalous interface diffusion effect is in this case? With q((p )A((j>o) — at x < — £d 
and q(4>o)A((f)o) = 1 at x > —£,d, the mass balance condition at the thin interface limit becomes 

V = -D L [d n u+ -Ud 2 s u], (10) 

where V is the interface velocity, u + is the diffusion field at the liquid side of the interface. The second term in the 
right hand side of Eq. (|l^) is the anomalous diffusion flux along the arc length of the interface. This flux originates 
from putting q(<po) A(<po) — 1 at — £d < x < within the solid side of the interfacial region, and the anomalous 
interface diffusion occurs within the half interface width ^d- With the choice 2£<j = Ax as explained before, this width 
corresponds to Ax/2. Whether or not the anomalous interface diffusion within the half mesh size affects much the 
dynamics of the interface is dependent on the computational systems and conditions. The effect can be significant 
when the curvature gradient inducing the interface diffusion is very large. Such situation may be met at the tip region 
during computations of dendritic growth with a large mesh size, which will be tested in the later part of this study. 
At the thin interface limit and under the conditions D$ <C Dl and q(^> )A(^o) = 1 at x > — following Ref. 
we find the velocity-dependent Gibbs-Thomson equation u = —d k — (3V with the kinetic coefficient 

where the constant a 2 is given by 

aiLx s(l)- 5 (-l) y/2f(M 

which appears to be identical with that in Ref. [UJtJ if the constant K therein is integrated partially, and also can be 
found in our previous study H by putting q((f>o)A(4>o) = 1 at the interfacial region. As pointed out in Ref. 0, Eq. 
( |i"T| ) and ( [i~2] ) remain unchanged in the symmetric case also where q((f>) = 1 and A(4>) = 1 in Eq. (||). 

III. COMPUTATION 

In present model, we adopted the parabolic potential /(</>) = |1 — </> 2 |/2 with cusps at <fi — ±1, the step-like function 
h(4>) defined by Eq. (||) and the discontinuous function q((j>) satisfying q{4>) = at <j> > (f>* (equivalently, at x < — £d) 
and q((f)A{(f)) = 1 at < cf>* . The combination of such functions in the model looks apparently unsuitable in view of 
the continuum spirit of PFMs. However the model works remarkably well in numerical computations to be shown. 

We simulated two cases. The first simulation was done on the dendritic growth in the symmetric case with q{4>) = 1 
(Ds = Dl) and A(<fi) = 1 (k = 1). We tested the effect of the interface width 2^ on the accuracy in tip velocity, as 
well as the performance of the parabolic potential. The second simulation was done on the one-sided case with Ds = 
and k < 1. We tested the convergence with do/W and the self-consistency in tip composition, and compared the tip 
velocity with that from Karma's model Q] . We solved the diffusion equation (||) and the well-known anisotropic form 
H of Eq. (0) by the typical explicit finite difference scheme. The four fold anisotropy in the interface energy was 
introduced by putting W(9) = W[l + £4 cos 40], where 9 is the angle between the direction normal to the interface 



and a>axis. The phase- field mobility was determined at the vanishing kinetic coefficient condition in Eq. (11); 



t{0) = a 2 XW(9) 2 /D L . (13) 

The functions and parameters commonly used in both simulations are as followings; f(<f>) = |1 — </> 2 |/2, g(4>) — 
(3/2) (<£ - 3 /3), h(4>) given by Eq. (§), °i = n / 4 > «2 = 0.2637 and 0.4150 for <j>* = 0.2 and 0.95, respectively, 
W = 1, Dl = 1 and the mesh size Ax = 0.4. Note that the interface width for the phase field then is given by 
2£p = kW = 7.85Ax, whereas the interface width 2^ = 2VFsin _1 </>* for diffusion field is given by l.OOAx and 6.27Ax 
for (j>* = 0.2 and 0.95, respectively. The initial undercooling (supersaturation) of the melt was uo — —0.55. The 
solidification started from the quarter circle solid seed of the radius tq — 40e?o at one corner of the system. The initial 
fields before solidification were u — 0, 4> — 1 at r < r and u — u — u Q exp[— e? (r — r )/(4Aa;)], (f> = — 1 at r > r . 

Two points must be explained for the numerical computations. The first is the stability problem related with the 
parabolic potential. At the edges x = ±£ p = ±irW/2 of the interface during computation, the phase-field values 
can be <f> > 1 or < — 1 and then oscillate around <j> — ±1 as computation proceeds, which is due to e?^/(±l) 7^ 
and the cusps at <j) = ±1 in the parabolic potential. Those oscillations propagate into the both the bulk solid and 
liquid region to make the interface unstable. This instability could be avoided by simply putting (j) = 1 if (j) > 1 and 
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<fi = — 1 if <j> < — 1 on all grids in every time steps. Even with such control, however, the interface became unstable 
when we took 4>* = 1 in Eq. @, that is h((p) — 4>, which is related with the anomalous release of latent heat (or 
solute redistribution) at the edges of the interface, following the artificial cutting of the phase-field values. For the 
computational stability, thus the adoption of the parabolic potential need to be combined with the localization of 
h{4>), as (f>* < 1 in Eq. (|J). The second point is the concentration oscillation in the solid of the one-sided system, 
which originates from the combination of the highly localized h((f>) and the discontinuous q(4>). After jumping from 
the interfacial state with — <f>* < <fi < <j>* to the solid state with <j> > <j>* at a given time step, the concentration of 
the grid is frozen. Thus the frozen concentration of the solid is dependent on the composition of the grid at the 
last time step remaining as the interfacial state. As the result, the frozen composition of the solid oscillates with 
solidification length and the oscillation amplitude becomes larger with decreasing number of grids in — ^ < x < £<j- 
At the dendrite tip of the one-sided system, the oscillation amplitude in most computations with 2^ = Ax was about 
10% of the composition shift due to Gibbs-Thomson effect. This composition oscillation could be removed by putting 
a small margin <5</> in the definition of g(</>); = 1 at 4> < (f>* + 5(f) and q(<fi) = otherwise. The minimum margin 
fi<fimin required to remove the oscillation corresponds to just the phase-field change during one time step in a grid 
with (f> = <fi* . In most computations of present study, <50 m in was less than 0.01. Therefore we adopted S(f> = 0.01 in all 
computations on the one-sided system, with which the oscillation amplitude in computations with 2£j = Ax could be 
suppressed within 0.1% of the composition shift due to Gibbs-Thomson effect. 

In the dendrite growth in undercooled pure melt (symmetric system) , the effect of the localization of latent heat 
release was tested. We took k = 1, e 4 = 0.05 and d /W =0.277 on 300 x 1200 grid system. In Fig. |l], the scaled tip 
velocity VcIq/Dl was plotted with the scaled time tD^/do for two different 4>* values in Eq. The thin line and 
dotted line are the results with <fi* — 0.2 and <f>* = 0.95 from this study, respectively, the thick line from the standard 
PFM, and the horizontal line is the steady tip velocity from the Green function method. For the standard PFM, we 
used the usual forms /(</>) = (1 - 2 ) 2 /4, g(cf>) =<j>- 20 3 /3 + 5 /5, h{4>) = 0, and then ai = 0.8839 from Eq. (§) and 
0,2 = 0.6267 from Eq. ([12]). As mentioned before, the case with 4>* — 1 could not be tested because the combination 
of the parabolic potential with the corresponding h(4>) — <f> makes the phase-field profile unstable. As can be seen in 
Fig. 0, two extreme interface widths of 2£d = Ax (0* = 0.2) and 2£d ~ 2£ p (0* = 0.95) for the diffusion field yield 
the similar tip velocity change with time. In particular, the tip velocity vs time for <jj* = 0.95 appeared to be nearly 
identical with that from the standard PFM h(<fi) = <f>. The final scaled tip velocities at steady state were 0.0172, 
0.0169, 0.168 for <f>* — 0.2, <f>* — 0.95 and the standard PFM, respectively, which are very close to the velocity 0.0170 
from the Green function method. Thus it turns out that the combination of the parabolic potential and the highly 
localized h(cj)) works as well as the standard PFM. 

Next we present the computation results for one-sided system. Fig. || shows the variation of scaled tip velocity 
VcIo/Dl with the scaled time tD^/d^, where we took k = 0.15 and 64 = 0.02. The thin line and the dotted line are 
the results for do/W = 0.277 on 300 x 1200 grid system and d /W = 0.554 on 600 x 2000 grid system, with the same 
4>* = 0.2. The dashed line is the result for do/W — 0.277 and <f>* = 0.95 on 300 x 1200 grid system. For comparison, 
the computational results obtained from the anti-trapping current model was included as the thick line in Fig. |^. 
Following the notations in this study, his model can be written as the phase-field equation (|l|) with /(</>) = (1 — 2 ) 2 /4, 
g(cf>) = <fi — 2(j> 3 /3 + cf> 5 /5, h(4>) = (j> and the diffusion equation 



The phase-field mobility at the vanishing kinetic coefficient condition is given by Eq. (|i3|), where a\ — 0.8839 and 

a 2 = 0.6267. As shown in Fig. the tip velocity variations with time remain unchanged for twice change in d /W 
when 2£d = Ax (or <f>* = 0.2) was adopted. If the surface diffusion could not be negligible, then its effect and so 
the tip velocity would be dependent on the grid size or do/W. This convergence thus indicates that the anomalous 
surface diffusion was negligible even at do/W = 0.277. The tip velocities for <jf = 0.2 are in close agreement with 
that (thick line in Fig. ||) from the anti-trapping model Q. We also measured the tip radius variation with the scaled 
time, following the method suggested in Ref. The tip radii for cf>* = 0.2 were within 5% error, compared with 
that computed from the anti-trapping model. These good agreement in tip velocity and tip radius may be seen as 
the evidence that the anomalous interface phenomena were effectively suppressed in both approaches. We increased 
the interface width for the diffusion field over six times from 2£d = Ax ((/)* — 0.2) to 2£<j = 6.27Ax (0* = 0.95), 
while keeping d /W — 0.277. The tip velocity then was decreased significantly as seen from the dashed line in Fig. |2j, 
whereas the tip radius was increased by 40%. Considering that the surface diffusion current decreases the curvature 
gradient of the dendrite tip, these must be the manifestation of the significant anomalous surface diffusion effect. Fig. 
|3| shows a steady-state composition profile across the interface along the growth axis of the dendrite in the one-sided 
alloy system with €4 = 0.02, k — 0.15, do/W — 0.277 and (jf = 0.2. The scaled solid composition for plane-front 
interface at equilibrium state is c e s /c e L — 0.15. Note that in Fig. ||tlie solute redistribution around the interface occurs 



-[uA{<f>) - ~h(<t>)} = • (1 - + ^=V ■[! + (!- *0«]j^r6 



7- 



(14) 
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within one grid width only due to the localization of h{<j>). The inserted magnified figure shows the comparison of 
the measured solid composition (filled circles) with prediction from Gibbs-Thomson equation (thick horizontal line), 
u = —d /p or cquivalently cs/c e L = k[l— (1 — k)d /p], where p is the computed tip radius. Excellent agreement within 
3 % error can be noticed between the computation results and the prediction. 

IV. SUMMARY 

In summary we minimized the interface diffuseness in the phase-field models by introducing the parabolic double- 
well potential and localizing the solute redistribution (or latent heat release) into a narrow region within a phase-field 
interface. In spite of the parabolic potential with cusps, highly localized solute redistribution and discontinuous diffu- 
sivity in this model, it works remarkably well in numerical computations. The computations on dendritic solidification 
of an one-sided system yield quantitatively the same results with Karma's anti-trapping model, indicating the anoma- 
lous interfacial effects can be effectively minimized. This approach can be easily extended to the multi-components 
or multi-phases system. This approach also is useful in suppressing the anomalous interaction between interfaces in 
the closely-spaced multi-particles systems, e.g. liquid-phase sintering. 
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FIG. 1. Effect of h((f)) width on the scaled tip velocity variations with the scaled time. The computational and material 
parameters were k = 1, Mo = —0.55, £4 = 0.05, W = 1, Dl = Ds = 1, Ax/W = 0.4 and the phase-field mobility r was 
determined at the vanishing kinetics condition. 

FIG. 2. Scaled tip velocity variations with the scaled time for one-sided alloy model with Ds = 0. The computational and 
material parameters were k = 0.15, uo = —0.55, £4 = 0.02, W = 1, Dl = 1, Ax/W = 0.4 and the phase-field mobility r was 
determined at the vanishing kinetics condition. 

FIG. 3. The composition profile across the interface along the growth axis of the dendrite and comparison of the measured tip 
composition with Gibbs-Thomson equation along the growth axis for one-sided alloy model with Ds = 0. The computational 
and material parameters were k = 0.15, u = -0.55, e 4 = 0.02, d /W = 0.277, W = 1, D L = 1, Ax/W = 0.4 and the 
phase-field mobility r was determined at the vanishing kinetics condition. 
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